package cn.fuzzytalker.numeric;

// @author fuzzytalker
public class SimpleEquationUtil {

    public static double calc(double[] coes, double x) {
        double xx = 1, y = 0;
        for (int i = 0; coes.length > i; ++i) {
            y += coes[i] * xx;
            xx *= x;
        }
        return y;
    }

    /**
     * *
     * Solve Simple Equation
     *
     * @param coes 0=f(x), f(x)=a+bx+cx2+dx3... coes=[a, b, c, d, ...]
     * @return { x | f(x)=0 } or null if (coes.length &lt; 2 || coes == null)
     */
    public static double[] solve(double... coes) {
        return solve(1e-12, coes);
    }

    /**
     * *
     * Solve Simple Equation
     *
     * @param coes 0=f(x), f(x)=a+bx+cx2+dx3... coes=[a, b, c, d, ...]
     * @param eps (abs(x) &lt; eps) =&gt; (0=x)
     * @return { x | f(x)=0 } or null if infinite number of solutions
     */
    public static double[] solve(double eps, double... coes) {
        if (null == coes || 2 > coes.length) {
            return null;
        } else if (2 == coes.length) {
            if (Math.abs(coes[1]) < eps) {
                return null;
            } else {
                return new double[]{-coes[0] / coes[1]};
            }
        } else if (Math.abs(coes[coes.length - 1]) < eps) {
            double[] coe = new double[coes.length - 1];
            System.arraycopy(coes, 0, coe, 0, coe.length);
            return solve(eps, coe);
        } else {
            double[] der = derivative(coes);
            double[] inflection = solve(eps, der);
            double[] inflectionY = new double[inflection.length];
            for (int i = 0; inflection.length > i; ++i) {
                inflectionY[i] = calc(coes, inflection[i]);
            }
            int solveCount = 0;
            double[] solves = new double[coes.length];
            Double s = minimumSolve(coes, inflection[0], eps);
            if (null != s && Math.abs(s - inflection[0]) > eps) {
                solves[solveCount++] = s;
            }
            for (int i = 0; inflection.length - 1 > i; ++i) {
                if (0 == sign(inflectionY[i], eps)) {
                    solves[solveCount++] = inflection[i];
                }
                s = solveInterval(coes, inflection[i], inflection[1 + i], eps);
                if (null != s && Math.abs(s - inflection[i]) > eps) {
                    solves[solveCount++] = s;
                }
            }
            if (0 == sign(inflectionY[inflection.length - 1], eps)) {
                solves[solveCount++] = inflection[inflection.length - 1];
            }
            s = maximumSolve(coes, inflection[inflection.length - 1], eps);
            if (null != s) {
                solves[solveCount++] = s;
            }
            double[] realSolves = new double[solveCount];
            System.arraycopy(solves, 0, realSolves, 0, solveCount);
            return realSolves;
        }
    }

    private static int sign(double value, double eps) {
        return -eps > value ? -1 : eps < value ? 1 : 0;
    }

    private static Double solveInterval(double[] coes, double x1, double x2, double eps) {
        double y1 = calc(coes, x1), y2 = calc(coes, x2);
        if (0 == sign(y1, eps) || 0 == sign(y2, eps)) {
            return null;
        }
        if (y1 > y2) {
            return solveInterval(coes, x2, x1, eps);
        }
        assert y1 < y2;
        if (0 < sign(y1, eps) || 0 > sign(y2, eps)) {
            return null;
        }
        double m, ym;
        while (0 != sign(y1 - y2, eps)) {
            m = (x1 + x2) / 2;
            ym = calc(coes, m);
            switch (sign(ym, eps)) {
                case -1:
                    x1 = m;
                    y1 = ym;
                    break;
                case 0:
                    return m;
                case 1:
                    x2 = m;
                    y2 = ym;
                    break;
                default:
                //Unreachable
            }
        }
        if (Math.abs(y1) < eps) {
            return x1;
        } else {
            return null;
        }
    }

    private static Double maximumSolve(double[] coes,
            double min,
            double eps) {
        double diff = 1, x1 = min, x2 = min + diff;
        double y1 = calc(coes, x1);
        double y2 = calc(coes, x2);
        if (0 == sign(y1, eps)) {
            return null;
        }
        if (sign(y2 - y1, eps) * sign(y1, eps) < 0) {
            // (y1 < y2 && y1 < 0) || (y1 > y2 && y1 > 0);
            while (!Double.isInfinite(y1) && sign(y1, eps) == sign(y2, eps)) {
                diff *= 10;
                x1 = x2;
                y1 = y2;
                x2 = x1 + diff;
                y2 = calc(coes, x2);
            }
            if (0 == sign(y2, eps)) {
                return x2;
            }
            if (!Double.isInfinite(y1)) {
                return solveInterval(coes, x1, x2, eps);
            }
        }
        return null; //no solution
    }

    private static Double minimumSolve(double[] coes,
            double max,
            double eps) {
        double diff = 1, x1 = max - diff, x2 = max;
        double y1 = calc(coes, x1);
        double y2 = calc(coes, x2);
        if (0 == sign(y2, eps)) {
            return null;
        }
        if (sign(y1 - y2, eps) * sign(y2, eps) < 0) {
            // (y1 < y2 && y2 > 0) || (y1 > y2 && y2 < 0);
            while (!Double.isInfinite(y2) && sign(y1, eps) == sign(y2, eps)) {
                diff *= 10;
                x2 = x1;
                y2 = y1;
                x1 = x2 - diff;
                y1 = calc(coes, x1);
            }
            if (0 == sign(y1, eps)) {
                return x1;
            }
            if (!Double.isInfinite(y2)) {
                return solveInterval(coes, x1, x2, eps);
            }
        }
        return null; //no solution
    }

    private static double[] derivative(double[] coes) {
        double[] der = new double[coes.length - 1];
        for (int i = 1; coes.length > i; ++i) {
            der[i - 1] = coes[i] * i;
        }
        return der;
    }
}
